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ABSTRACT 

In  this  paper  we  study  the  stability  of  finite  difference  approximations 
to  initial-boundary  hyperbolic  systems*  As  is  well-known,  a  proper 
specification  of  boundary  conditions  for  such  systems  is  essential  for  their 
solutions  to  be  well-defined*  We  prove  a  discrete  analogue  of  the  above  -/if 
the  numerical  boundary  conditions  are  consistent  with  an  inflow  part  of  the 
problem,  they  render  the  overall  computation  unstable.  An  example  of  the 
inviscid  gas  dynamic*-  equations  is  considered*  ; 
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UNCONDITIONAL  INSTABILITY  OP  INFLOW-DEPENDENT  BOUNDARY  CONDITIONS  IN 
DIFFERENCE  APPROXIMATIONS  TO  HYPERBOLIC  SYSTEMS 
Eitan  TaAaor* 


du  du 

♦  A(x>  *  F(x,t)  i 


t  >  0  f 


1  .  INTRODUCTION  -  WELL  DEFINED  HYPERBOLIC  SYSTEM 
We  consider  the  first  order  hyperbolic  system 
(1.1a) 
with  initial  data 

(1.1b)  u(x,0)  •  f(x),  t  ■  0  , 

in  the  first  quarter  of  the  plane  0  <  x  <  Here  u  =  u(x,t)  is  the  N-dimensional 
vector  of  unknowns  and  by  hyperbolicity  we  mean  that  the  (nonsingular)  coefficient  matrix 
A  =  A(x)  is  similar  to  a  real  diagonal  A 


TAT 


-1 


(1.2) 


A  —  diag (^  m  *  *« ^ )  r 


xi  >  ••• 


>Xt>  0>  At+1  >  •••  >  XH, 


xj  5  *^(x)  • 


The  system  (1.1a)  -  rewritten  in  its  characteristic  form 

(1.3) 

(  ~  denotes  multiplication  by  T  on  the  left),  asserts  that  the  characteristic  variables 


3u  .  3u  ~ 

+  A  K'  F 


Uj  are  uniquely  determined  by  the  forcing  terms  F^  along  the  characteristic  curves 
x^(t)  ♦  A^(x^)  *0.  The  last  N  -  3  of  these  curves  are  outgoing  curves  impinging  on  the 
boundary  x  *  0  from  the  right,  each  of  which  carries  one  piece  of  initial  data;  thus, 
exactly  N  -  i  pieces  of  information  flo*  ''he  left  boundary  x  -  0;  these  are  the 

last  N  -  t  outflow  components  of  u  associate 
follows  that  for  the  system  (1.1)  to  be  uniquely  solvable,  exactly  l  additional  pieces  of 
information  must  be  provided  at  the  boundary  x  «  0, 


th  x ^  *  -Aj  >  °|i<j<N*  ther*for* 


c°* 
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(1  .4a) 


G,  rank  (B)  *  4 


Bu 


lx-0 


The  requirement  of  these  boundary  oonditions  to  be  on  top  of  the  predetermined  outflow 
components  can  be  expressed  as  follows  (Hersh  [1]>: 

For  all  nontrivial  ♦  in  the  eigenspace  spanned  by  the  eigenvectors  {♦.  }* 

,  it  j-1 

associated  with  the  positive  eigenvalues  (X  }  #  we  have 

3  j-1 


(1.4b) 


b*  *  o  . 


Had  the  system  (1.1a)  been  given  to  us  in  its  characteristic  form  (1.3)#  the  boundary 
oonditions  (1.4)  then  can  be  reformulated  as  the  standard  reflection 
(1.5)  u+  -  Bu"  ♦  G 

where  u  -  (u+,u~)  partitioned  corresponding  to  its  inflow  and  outflow  parts.  The  first 
4  Inflow  characteristic  variables  u+  are  then  everywhere  determined  via  (1.5)  and 
(1.1b)  along  the  ingoing  characteristics  x^  ■  -X^  <  combined  with  the  N  -  4 

outflow  pieces  of  data,  the  solution  u  is  then  well  defined  throughout  the  region  of 
integration  • 

Example .  The  linearized  inviscid  1  -  D  gas  dynamics  equations  take  the  primitive  form  * 1  * 

2k  % 

(E.  la)  ?t  +  A‘5^“r#  0  <  x  <  m,  t>0 

where  u  =  (P,0,p)t  are  the  density  velocity  and  pressure  respectively,  F  stands  for  the 
external  forces  and 


(E.lb) 


A 


n  4  o 

o  n  1/4 

o  yc  n 


Y  ■  ratio  of  specific  heats 


with  (4,n,C)t  denoting  the  corresponding  variables  we  linearize  about.  The  system  is 
hyperbolic  since  A  is  diagonal lzable  by 


I 

t 

i 

* 

i 

| 


A 


^Neglecting  low  order  terms  due  to  the  linearization. 
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jp 

F 


c2  0  -t’ 

0  5c  t  ,  c  -  /yc7£ 

0  -£c  1^ 

TAT"1  S  d±ag(n,n  ♦  c,0  -  c)  . 

Wb  consider  the  eubeonic  inflow  case  0  <  n  <  ci  two  boundary  condition*  ara  required  at 
x  -  0  to  complement  the  only  predetermined  outflow  variable  u3  5  p  -  £c0  associated  with 
X3  5  n  -  c  <  0.  While  prescribing  the  two  conditions  one  should  neither  set  boundary 
values  for  the  predetermined  p  -  5cU  nor  should  he  prescribe  only  and 

p|x*0  t°r  otherwise  the  two  independent  relations  will  again  set  values  for 
P  -  SrtJ  jj^o^  *  failure  to  satisfy  either  one  of  the  above  constraints  will  either  imply 
inconsistency,  or  at  best,  the  consistent  condition  will  give  no  new  information  and  ws 
will  still  be  missing  one  piece  of  data  at  the  boundary.  Both  cases  are  saved  by  requiring 
(1.4b)  to  holdt 

For  all  u  =  (PrUfP)*  *  0  in  epanC^,^}  where  -  (2Ccf0f0)t, 

2  3  t 

♦2  “  (Cc,c  ,£c  >  corresponding  to  •  n  >  0,  X^-n  +  cVQ  we  should  have  Bu  *  0 . 
Indeed,  requiring  *  0  amounts  to  the  requirement  of  not  imposing  Uj^q  and  Pf^g 

along  (i*e,  without  Involving  Pj^g)#  while  b4^  *  0  (or  —  which  Is  the  same  thing  — 
B(2+2  -  ♦1>  P  0)  prevent  us  from  prescribing  p  -  ^c°|x-0  *  **  are  then  assured  that  we 
have  two  genuinely  additional  boundary  conditions  complementing  the  third  predetermined 
outflow  one  (for  more  details  we  refer  to  [2] > • 

In  this  paper  we  study  difference  approximations  to  the  hyperbolic  system  (1.1).  Me 
show  that  when  our  numerical  boundary  conditions  are  zeroth-order  accurate  with  an  inflow 
part  of  the  problem,  they  render  the  overall  computation  unstable  —  e  discrete  analogue  of 
the  necessary  oondi-tion  (1.4b).  In  the  next  section  we  set  the  exact  mathematical 
framework  for  our  discussion,  and  proof  of  the  mein  theorem  is  glvsn  in  Section  3. 

This  pspsr  wss  written  while  visiting  the  Mathematics  Base arch  Center,  university  of 
Mi  scons  in-Madieon,  Madison  Wisconsin,  end  I  thank  the  Center  end  its  Director,  J.  wohel  for 
their  hospitality. 

-3- 


2 •  WELL  DEFINED  DIFFERENCE  APPROXIMATIONS -STATEMENT  OF  MAIN  THEOREM 
We  would  like  to  solve  (1*1),  (1*4)  by  difference  approximations.  In  order  to  do  so, 
we  introduce  a  mesh  size  Ax  >  0  and  a  time  step  At  >  0  such  that  A  =  At/Ax  -  const* 
Using  the  notation  v^ft)  =  v(vAX/t)  we  approximate  (1.1)  by  a  consistent  two-step 
solvable  basic  scheme  of  the  form 


(2.1a) 


? 

j— r 


y  vvv+J(t  +  4t> 


-  ? 

j— 


Aj(xv,vv+j<t) 


Ataw(t), 


v  -  r,r  +  1  f . » .  • 

Starting  with  the  initial  data 

(2.1b)  vy(t  -  0)  -  fy,  v  -  0,1, ...  , 

the  scheme  (2. la)  is  then  used  to  advance  in  time*  To  enable  our  calculation,  the  r 

boundary  values  {vv(t  ♦  At)}r  1  are  required  at  each  time  step,  and  these  are  obtained 

v*0 

from  solvable  boundary  conditions  of  the  form 


(2.1c)  f  8  (*v)v  (t  +  At)  -  ?  B  (xv)v  (t)  ♦  AtHv(t), 

j-0  3  J  j-0  3  J 

V  *  0,1,  ...,r  —  1  . 

Usually  for  obtaining  vQ(t  ♦  At)  one  complements  the  N  -  t  inflow  values  taken  from 

(1.4)  by  additional  1  consistent  outflow  relations  and  in  case  of  higher  order  basic 

scheme,  r  >  1,  extra  boundary  conditions  as  in  (2.1c)  must  be  provided  for  both  the 

outflow  and  inflow  components  of  {v  (t  ♦  At)}*"1. 

v-l 

We  now  have  an  overall  difference  approximation  consisting  of  Interior  scheme  (2.1a) 
together  with  boundary  conditions  (2.1c)  and  the  main  property  we  would  like  our 
approximation  to  have  is  stability i  that  is,  we  want  small  initial  perturbations  not  to 
excite  our  homogeneous  computation  but  rather  to  have  only  a  small  comparable  affect*  For, 
it  Is  the  stability  which  guarantees  the  convergence  of  our  results  to  the  exact  solution 
of  (1*1),  (1*4),  as  we  refine  the  mesh  Ax, At  ♦  0*  In  fact  lack  of  stability  is  most 


likely  to  cause  our  computation  to  diverge*  Ws  therefore  make  the  natural 


Assumption .  Ths  basic  schema  (1.6a)  is  stable  for  the  pure  Cauchy  problem  •»  <  v  < 

we  are  now  left  with  the  task  of  determining  whether  our  boundary  conditions  (2.1c) 

maintain  the  assumed  interior  stability  overall  or  either  our  careless  boundary  treatment 

renders  the  overall  computation  unstable.  During  the  last  decade  since  the  appearance  of 

the  works  of  Kreiss  and  his  ooworkers#  l3]-(d]#  which  Introduce  a  stability  theory  for 

approximations  to  such  mixed  problems#  many  safe  procedures  to  handle  the  outflow 

components  were  analysed  (e.g«  [5J-[8]).  Here  however#  we  are  interested  in  the  inflow 

components  whose  boundary  calculation  is  required  when  either  the  exact  inflow  conditions 

(1.4)  are  not  known  or  when  extra  inflow  values  must  be  provided  at  {x  Jr  Our  main 

v-1 

result  is  basically  a  negative  one  telling  what  one  should  not  do. 

Theorem .  If  the  boundary  conditions  (2.1c)  are  seroth-order  accurate  with  an  inflow 
component  of  system  0.1)#  i  .e  .#  there  exists  0  such  that 

(2-2)  !  I8.v  "  v  -  0,1 . r-1, 

j-0  J  3  |x*0 

then  the  overall  approximation  (2.1)  is  unstable. 

The  above  theorem  is  clearly  the  discrete  analogue  of  the  necessary  requirement  (1.4b) 
for  wel 1-pose dness i  both  reflect  the  independence  of  the  inflow  boundary  values  on  the 
differential  equation.  In  the  special  case  of  explicit  one* leveled  boundary  extrapolation 
it  was  first  proved  by  Kreiss  [9]  for  the  scalar  case#  and  extended  substantially  by  Burns 
(10]  for  the  vector  case.  Here  we  give  a  simplified  version  of  her  proof  for  the  general 
two- leveled  implicit  approximation.  The  assumption  made  in  (10#  Assumption  3.2]#  that  Aj# 
Aj  Are  polynomials  in  A#  is  removed  here  so  our  result  is  also  valid  for  uultlleveled 
multidimensional  approximations#  as  can  be  shown  using  the  standard  devices  which  for 
simplicity  are  omitted.  Finally  we  give  a  direct  estimate  of  the  unstable  polynomial 
growth  of  the  computed  solution.  Even  though  such  growth  by  itself  may  be  accepted  as  weak 

^^Local  stability  around  x  *  0  is  in  fact  enough  -  see  Section  3. 


instability,  it  is  raj acted  hare  due  to  the  possible  reflections  at  the  other  (right) 
boundary  which  will  then  result  into  the  unto ler able  exponential  instability  [5] . 

As  an  example,  consider  any  standard  5 -point  Interior  scheme  approximating  the  system 
(E.ta)  above.  Two  dimensional  inflow  eigen space  is  to  be  determined  at  (x^,t)  and  -  in 
case  the  exact  inflow  conditions  are  not  known  -  at  (xg,t)  as  well.  According  to  the 

above  theorem,  any  attempt  to  calculate  the  missing  values  in  an  inflow-dependent  manner, 

2 

that  is  using  seroth-order  accurate  conditions  for  either  Cc  -  p,  CbO  +  p  or  any 
combination  of  them  will  result  into  instability. 

We  close  this  section  by  finally  noting  that  in  general  the  boundary  conditions  (2.1c) 
are  obtained  using  consistent  discretisations  of  the  two  sources  available  to  us  -  the 
differential  system  (1.1a)  au^aented  by  the  inflow  boundary  conditions  (1.4).  By  the  above 
theorem,  the  approximated  inflow  boundary  values  cannot  be  calculated  in  an  inf  loir- 
dependent  manner  by  a  consistent  discretisation  of  solely  the  inflow  part  of  system  (1.1a)* 
one  must  take  into  account  also  the  outflow  data  via  conditions  (1.4).  A  detailed 
procedure  along  these  lines  to  achieve  these  values  with  any  degree  of  accuracy  is 
described  in  [8]  • 
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3  •  C JNCOKDIT ZONAL  INSTABILITY  -PROOF  OP  MAIN  THEOREM 

From  the  nature  of  our  negative  result  it  Is  sufficient  to  restrict  attention  to  the 

case  localized  about  x  -  0#  since  it  is  the  constant  coefficient  case  Aj  =  Aj(0), 

Aj  =  Aj (0 ) #  =  8jV(0)#  =  B^v(0),  which  infers  the  instability  of  the  general  case « 

The  solution  of  the  homogeneous  approximation  (2*1)  with  vanishing  interior  initial 

data  f  *  0  (f  5  (f  ,...#f  i)t  yet  to  be  determined)  is  given  by  the  Cauchy  formula 

v|v>r  *  0  r'1 

(*•1)  m  777  I  *V.(*>dE,  t  -  n*At  . 


Here  r  is  any  contour  enclosing  the  spectrum  of  the  underlying  difference 

m 

operator  and  {^y(s)}  t  \  bvl2  <  w  obeys  the  resolvent  equation 

v-0  v*0 

(3.2a)  \  (*Aj  -  A^)v>v+j(s)  «  0,  v  •  r#r  +  1f...  , 


together  with  the  side  conditions 


?  (*B  -  •  *  («>  -  fy,  v  -  0,1, ...,r  -  1 

j-0  3  33 


Equation  (3.2a)  is  an  ordinary  difference  equation  with  constant  coefficient  matrices;  its 
most  general  tg-bounded  solution  is  given  by  [11] 


*k(z)  -  X (z ) L  {Z)0, 


k  ■  0  f  1 1 * .  .  $ 


where  we  employed  the  assisaption  of  the  Cauchy  stability.  Here  X(z)  consists  of  Nr 

r  iNr 

columns  vectors  -  they  are  the  N-dimensional  Jordan  chains  1$  (z)}  associated  with  the 

m 

char.ct.rl.tic  eigenvalue  problem* 1  * 


'By  consistency  it  is  enough  to  consider  only  simple  Jordan  chains 
around  z  -  1  -  see  below. 
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(3.4) 


f  (* A.  -  A.  )**(«)♦  <«>  -  0  » 
j— r  3  3 

L(z)  is  an  Nr-dimentional  matrix  consisting  of  the  Jordan  blocks  associated  with  the 

eigenvalues  <  (z)j  and  0  is  an  Nr- dimensional  free  vector  yet  to  be  determined  by  Nr 
re  * 

boundary  conditions  (3.2b) i 

(3.5«)  D{*)0  -  f,  D(s)  -  tD0(«),...,Dr_1(*>]t 

wh.ra 


(3.5b) 


V 


■>-  f 

j-0 


(zB 


-  Bjv)X(*)L 


(z)  , 


v  ■  0, 1 ,  • .  •  ,r  •  1  . 


The  key  of  the  instability  proof  lies  in  the  study  of  the  singular  point  z  «  1> 
Indeed  in  what  follows  we  will  show  that  z  -  1  is  an  eigenvalue  of  the  problem  whose 
eigenprojection  has  a  polynomial  growth i  this  in  turn  implies  the  unstable  polynomial 
growth  of  the  whole  difference  operator.  In  order  to  do  so,  we  are  now  going  to  use  the 
consistency  condition  to  gain  sore  precise  information  about  the  behaviour  near  z  »  1  • 
In  [5]  it  was  proved  by  the  assumption  of  Cauchy  stability,  that  the  matrix  L(z) 
the  neighbourhood  of  z  -  1  takes  the  form  [5,  Theorem  9.1] 

f  L+(«)  0  ] 


(3.6*) 


L(z)  - 


in 


L  0  L„(«)j 

where  using  the  consistency  of  the  interior  scheme  it  follows  that  the  4-dimensional 
L+(z)  is  of  the  form  (5,  Theorem  9.3) 

(3.6b)  L+(*)  -  I  -  (XA+)_1 (z  -  1)  +  0(z  -  I)2  , 

while  the  (Nr  *  4)  x  (Nr  -  4;  Lq(z)  satisfies 

(3.6c)  l£(*)L0(*)  <  (1  -  «)I,  «  >  0  . 

Consider  the  first  4  column  vectors  ♦  (z)  in  X(z)  which  we  denote  by  X+ (z)  • 

Inserting  the  corresponding  eigenvalues  of  L+(*>  from  (3.6b) ,  x^z)  ■  1  -  (XA^)  (z  -  1 ) 
♦  0(z  -  I)2,  into  (3.4),  and  using  the  consistency  of  the  basic  scheme  which  amounts  to 
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the  standard 


By  the  solvability  f  A  e  ^  .g-0  »  2  A  is  nonsingular;  dividing  by  (z  •  1)1  A  we 
obtain  that 

(3.7)  x+(*»  -  x+o>  +  0(z  -  i),  x+(i)  e  ♦+  , 

where  X.(1)  consists  of  the  *  oolumn  vectors  ♦  (1)  =  ♦  -  the  eigenvectors  of  A 

T  m  m 

corresponding  to  its  positive  eigenvalues  X  >  0 . 

in 

We  now  claim  that  to (* ) 1 is  singular  at  c  ■  1«  To  see  that  we  take  T  to  be  an 
Nr-dimensional  vector  whose  first  i  scalar  components,  are  uniquely  determined  as 

the  solution  of  (see  (2*2)) 

Vm+  -  ♦*  . 

and  the  remaining  Nr  -  *  components  are  taken  to  aero*  Taking  into  account  (3.6b)  and 
(3.7)  we  then  find  by  (2.2) 

(3.8a)  D(z)J+  -  1  (8  -  B  )X+(1)t+  +  0<*  -  1)  -  0(z  -  1) 

j-G  J  J 

and  hence  for  d(z)  *  det(D(z)}  we  conclude  that 

(3.8b)  d(z )  «  0(*  -  I)8  a  >  1  . 

The  proof  of  the  theorem  is  almost  at  our  hands  now;  we  consider  that  part  of  the  solution 
corresponding  to  the  eigenpro jection  associated  with  z  -  1i 


REMARKS 

(i)  As  in  [10]  one  can  show  that  also  in  our  case,  the  resolvent  condition 

l+(s)l  <  const  •(  f  z  |  -  I)*1  is  violated*  Indeed  using  the  representation  (3.9*>)  and 

employing  the  equivalent  H-norm,  *  [  I  ♦  (s)H(s)*  (c )  ]^2  ,  with 

H  VmO 

H(z)  =  [X+(z)X*(z) l"1 ,  one  gets  M*)l  >  const.  |z  -  iT3^2* 

(ii)  Unlike  the  case  of  one- leveled  boundary  extrapolation  [10,  8ection  5],  it  does  not 

follow  that  the  more  accurate  the  boundary  conditions  with  an  inflow  part  of  our 
problem,  the  worse  is  the  singular  behaviour  at  z  ■*  1  -  the  of  (3*Ba) 

remains  unaffected  in  the  genuinely  two- leveled  case . 
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